
;;==============================================================
;; Goal: Grid out T_A "detection" threshold & multiple-peak
;;   primary:secondary ratio to determine best values in terms
;;   of matching % with HCO+ VLSR & number of sources.

COMPILE_OPT IDL2, LOGICAL_PREDICATE
COMMON OMNI_CONFIG, conf, mw, local, dpdfs, ancil, fmt, conffile
COMMON GRS_BLOCK, grs, grsexist, v_std

goto,plot_only

conf = omni_load_conf()
s = omni_read_cat(conf.cat,ncat)
v_std = findgen(conf.nvbin)*conf.deltav + conf.vstart

IF ~exist(ancil) THEN $
   ancil = omni_read_conffile('./conffiles/ancillary.conf')

n_ta  = 81
n_rat = 80
dx = 1.d/n_rat

ta  = 10.^(dindgen(n_ta)/double(n_rat)*2.d - 2.d)       ; 0.01K -> 1K in LOG
rat = (dindgen(n_rat)+1.d)/double(n_rat)*2.d + 1.d + dx ; 1+dx -> 3

match = dblarr(n_elements(ta),n_elements(rat))
num   = dblarr(n_elements(ta),n_elements(rat))
count = bytarr(n_elements(ta),n_elements(rat))
fmult = dblarr(n_elements(ta),n_elements(rat))

;; Restore the velocity structure
restore,'./local/'+conf.survey+'_velocities.sav',/ver
IF n_elements(grs) EQ 0 THEN $
   restore,'./local/'+conf.survey+'_grs_spectra_r46.sav',/ver

start_t = systime(1)
openw,lun,'./grs_paper/tr_runtime.txt',/get_lun

randi = permute(n_elements(ta))
randj = permute(n_elements(rat))

;; Find where HCO+ velocities exist
ih = where(v.mol[0].nspec NE 0, nh)

FOR rr=0L,n_elements(ta)*n_elements(rat)-1 DO BEGIN
;; FOR ri=0,n_elements(ta)-1 DO BEGIN
;;    FOR rj=0,n_elements(rat)-1 DO BEGIN
   ri = rr MOD n_elements(ta)
   rj = rr MOD n_elements(rat)
   
   ii = randi[ri]
   jj = randj[rj]
   
   ancil.grs_ta = ta[ii]
   ancil.grs_ratio = rat[jj]
   
   print,'----------------------------------------------------------'
   print,'Indices: ',ii,jj
   print,'Grid points: ',ancil.grs_ta,ancil.grs_ratio
   
   ;; Compute GRS spectra for these value of TA & RATIO
   message,'Computing GRS spectra...',/inf
   
   ;; First, clear out existing data in v.grs
   v.grs.flag    = 0b
   v.grs.vlsr    = 0.d
   v.grs.tpk     = 0.d
   v.grs.lw      = 0.d
   v.grs.tonspec = 0.d
   
   FOR ci=0,ncat-1 DO BEGIN
      
      ;;***************************************************************
      ;; Add info from the GRS 13CO sepectra into the structure
      ic = where(grs.cnum EQ v[ci].cnum, nc)
      
      IF nc THEN BEGIN          ; Skip if can't find in GRS
         spectrum = grs[ic].spectrum
         IF total(spectrum) EQ 0. THEN CONTINUE ; If no GRS to begin with
         omni_grs_maskspec, spectrum, grsflag

         void = Check_Math()

         v[ci].grs.flag    = grsflag
         IF grsflag EQ 0 || grsflag EQ 3 THEN CONTINUE ; If flagged bad
         yf = mpfitpeak(v_std,spectrum,A,NTERM=3)
         v[ci].grs.vlsr    = A[1]
         v[ci].grs.lw      = A[2]
         v[ci].grs.tpk     = max(spectrum,grsj)
         v[ci].grs.tonspec = grs[ic].onspec[grsj]
         
      ENDIF
      ;;*************************************************************** 
   ENDFOR
   
   
   
   
   ;; Find where GRS spectrum exists
   ig = where(v.grs.flag NE 0 AND v.grs.flag NE 3,ng)
   print,m4_stat(v.grs.flag)
   print,'NG: ',ng
   print,'NH: ',nh
   
   ;;============================
   ;; Compute matching stats!
   
   ;; Get list of GRS-HCOP compatabilities
   iginh = where_array(ig,ih,ngh)
   ihing = where_array(ih,ig,nhg)
   print,'NGH, NHG: ',ngh,nhg
   ihp = ih[iginh]
   igp = ig[ihing]
   nmult = double(total(v[ihp].grs.flag EQ 2))
   
   ikms = where(abs(v[ihp].vlsr-v[ihp].grs.vlsr) LE 5., nkms)
   match[ii,jj] = double(nkms)/n_elements(ihp)*100.d
   num[ii,jj] = ngh
   fmult[ii,jj] = nmult/ngh*100.d
   print,'Stats: ',match[ii,jj],num[ii,jj],fmult[ii,jj]
   count[ii,jj] = 1b
   print,'Grid points completed: ',long(total(count))
   ;;============================
   
   
   ;; Figure time elapsed and ETD...
   elapsed = systime(1) - start_t
   done = double(rr+1) / double(n_elements(ta)*n_elements(rat))
   print,'ETD: ',systime(0,start_t + elapsed/done),$
         '  ['+string(rr+1,n_elements(ta)*n_elements(rat),$
                      format="(I0,'/',I0,']')")
   
   printf,lun,elapsed,elapsed/done
   flush,lun
ENDFOR

close,lun
free_lun,lun

save,match,num,fmult,count,n_ta,n_rat,$
     filename='./grs_paper/ta_rat_planes.sav',/ver


;;**************************************************************************
;;==========================================================================
;;**************************************************************************
plot_only:

IF ~exist(ancil) THEN $
   ancil = omni_read_conffile('./conffiles/ancillary.conf')

restore,'./grs_paper/ta_rat_planes.sav',/ver

;; Filter the arrays for easier plotting
match = median(match,5)
num = median(num,5)

;; Compute useful things for the plotting
dx = 2./n_ta
dy = 2./n_rat

help,match
ta_p  = 10.^(dindgen(n_ta)/double(n_rat)*2.d - 2.d) ; 0.01K -> 1K in LOG
rat_p = (dindgen(n_rat)+1.d)/double(n_rat)*2.d + 1.d + 1.d/n_rat ; 1+dx -> 3

ita  = where( abs(ta_p  - ancil.grs_ta)  LE dx/2., nita)
irat = where( abs(rat_p - ancil.grs_ratio) LE dy/2., nirat)
print,'Prints:'
print,nita,nirat
print,ancil.grs_ta,ancil.grs_ratio
print,ta_p[ita]
print,rat_p[irat]
print,match[ita,irat[0]:irat[-1]]
print,num[ita,irat[0]:irat[-1]]
print,'Prints'

;;====================================================================
;; PLOT!!1!

myps,'./grs_paper/ta_rat_optimization.eps',xsize=3.1,ysize=2.8,thick=1.5,$
     ct=34
!x.omargin = [-0.75,+2.5]
!y.omargin = [-0.5,-1.5]
multiplot,[1,1]

csf = 0.55
;; multiplot_xm,[2,1],mpcharsize=1.0,xgap=0.02,/hardleft,/hardright
cgLoadct,34

matlevel = [90,92,94,95,96,97,98,100]

;;====================================================================
;; Panel LEFT

pr = set_plot_range(num)
plotimage,num,range=pr,$;title='Number of GRS - HCO!u+!n(3-2) objects',$
          xtit='Thresh',ytit='Ratio',xst=4,yst=4,min_dpi=300

dx = 1./n_ta
dy = 2./n_rat
cgAxis,xaxis=0,/xst,xr=10.^[-2-dx/2.,dx/2.],/save,/xlog,$
       xtit='T!dA,thresh!n  [K]',charsize=1.0*csf
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,/yst,yr=[1+1.d/n_rat-dy/2.,3+1.d/n_rat+dy/2.],/save,$
       ytit='R!dT!n',charsize=1.0*csf
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

print,'Ranges? '
print,'X: ',[10.^[-2-dx/2.,dx/2.],minmax(ta_p)]
print,'Y: ',[[1+1.d/n_rat-dy/2.,3+1.d/n_rat+dy/2.],minmax(rat_p)]

;; Number of objects in the overlap
nobjs = [900,1000,1100,1200,1300,1400,1500]
cgContour,num,ta_p,rat_p,levels=nobjs,/over,color='blk7',c_charsize=1.0*csf,c_linestyle=5
;; cgContour,num,ta_p,rat_p,levels=[1035,1040],/over,color='deep_pink',c_charsize=1.0*csf,c_linestyle=5

;; Multiplicity Fraction
lvs = indgen(20)*10.
cgCOntour,fmult,ta_p,rat_p,c_annot=string(lvs,format="('   ',I0,'%   ')"),$
          color='black',/over,levels=lvs,c_charsize=0.9*csf

;; Matching percentage
cgContour,match,ta_p,rat_p,levels=matlevel,/over,c_charsize=1.0*csf,$
          c_color='wt1',c_annot=string(matlevel,format="(' ',I0,'% ')")




;; vline,0.05,color='wt1'

cgPlots,psym=14,color='blk5',symsize=1,0.30,1.7
cgPlots,psym=4,color='black',symsize=1,0.30,1.7
cgPlots,psym=14,color='blk5',symsize=1,0.06,2.05
cgPlots,psym=4,color='black',symsize=1,0.06,2.05
cgPlots,psym=14,color='white',symsize=1,ancil.grs_ta,ancil.grs_ratio
cgPlots,psym=4,color='black',symsize=1,ancil.grs_ta,ancil.grs_ratio

cgColorbar,/vert,/right,position=[0.87,!y.window[0],0.90,!y.window[1]],$
           range=pr,charsize=0.9*csf,title='Number of GRS - HCO!u+!n(3-2) objects'




cgAxis,yaxis=0,/yst,ytickformat='blank_axis'
 
;; multiplot,/doyaxis

;; ;;====================================================================
;; ;; Panel MIDDLE

;; print,'MATCH: ',m4_stat(match)

;; pr = [88,97]
;; plotimage,match,range=pr,title='Matching Percentage',$
;;           xtit='Thresh',ytit='Ratio',xst=4,yst=4,/preserve,min_dpi=300

;; cgContour,match,levels=matlevel,color='black',/over

;; cgAxis,xaxis=0,/xst,xr=[0+dx/2.,1+dx/2.],/save,$
;;        xtit='T!dA!n Threshold  [K]',charsize=1.0
;; cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
;; cgAxis,yaxis=0,/yst,yr=[1-dy/2.,3+dy/2.],/save,$
;;        ytit='Primary:Secondary Ratio',charsize=1.0
;; cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

;; cgPlots,psym=14,color='brown',symsize=2,0.3,1.5
;; cgPlots,psym=4,color='black',symsize=2,0.3,1.5
;; cgPlots,psym=14,color='blk5',symsize=2,0.2,1.5
;; cgPlots,psym=4,color='black',symsize=2,0.2,1.5
;; cgPlots,psym=14,color='bisque',symsize=2,0.1,1.55
;; cgPlots,psym=4,color='black',symsize=2,0.1,1.55
;; cgPlots,psym=14,color='blu3',symsize=2,0.05,1.60
;; cgPlots,psym=4,color='black',symsize=2,0.05,1.60

;; cgColorbar,/vert,position=[0.02,!y.window[0],0.035,!y.window[1]],$
;;            range=pr,charsize=0.5

;; multiplot,/doyaxis

;;====================================================================
;; Panel RIGHT

;; plotimage,fmult,range=[0,100],title='Fraction with Multiple Peaks',$
;;           xst=4,yst=4,/preserve,min_dpi=300
;; cgContour,match,levels=matlevel,color='wt1',/over,$
;;           c_annot=string(matlevel,format="(' ',I0,'% ')")
;; lvs = indgen(20)*10.
;; cgCOntour,fmult,c_annot=string(lvs,format="('   ',I0,'%   ')"),$
;;           color='brown',/over,levels=lvs

;; dx = 1./n_ta
;; dy = 2./n_rat
;; cgAxis,xaxis=0,/xst,xr=10.^[-2-dx/2.,dx/2.],/save,/xlog,$
;;        xtit='T!dA,thresh!n  [K]',charsize=1.0
;; cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
;; cgAxis,yaxis=0,/yst,yr=[1-dy/2.,3+dy/2.],/save,$
;;        ytit='R!dT!n',charsize=1.0
;; cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

;; ;; vline,0.05,color='wt1'
;; ;; cgPlots,psym=14,color='brown',symsize=2,0.3,1.5
;; ;; cgPlots,psym=4,color='black',symsize=2,0.3,1.5
;; ;; cgPlots,psym=14,color='blk5',symsize=2,0.2,1.5
;; ;; cgPlots,psym=4,color='black',symsize=2,0.2,1.5
;; ;; cgPlots,psym=14,color='bisque',symsize=2,0.1,1.55
;; ;; cgPlots,psym=4,color='black',symsize=2,0.1,1.55
;; cgPlots,psym=14,color='bisque',symsize=2,ancil.grs_ta,ancil.grs_ratio
;; cgPlots,psym=4,color='black',symsize=2,ancil.grs_ta,ancil.grs_ratio




myps,/done,/mp


;;================================================
;; Do some analysis...
print,'--------------------------------------------------'
ta  = 10.^(dindgen(n_ta)/double(n_rat)*2.d - 2.d)       ; 0.01K -> 1K in LOG
rat = (dindgen(n_rat)+1.d)/double(n_rat)*2.d + 1.d + dx ; 1+dx -> 3

tt = cmreplicate(ta,n_rat)
rr = transpose(cmreplicate(rat,n_ta))

help,match,tt,rr

ind = where( abs(match-95.) LE 0.1, ni)
nm = num[ind]

print,max(nm,jj)
print,tt[ind[jj]],rr[ind[jj]]


print,m4_stat(num[ind],tt[ind],rr[ind])

print,'The chosen one:'
ind2 = where( abs(alog10(tt) - alog10(ancil.grs_ta)) LE 0.007 AND abs(rr-ancil.grs_ratio) LE 0.03, ni2)
print,ni2

print,num[ind2],match[ind2],tt[ind2],rr[ind2]


print,'Off point 1:'
ind3 = where( abs(alog10(tt) - alog10(0.30)) LE 0.007 AND abs(rr-1.70) LE 0.03, ni3)
print,ni3

print,num[ind3],match[ind3],tt[ind3],rr[ind3]

print,'Off point 2:'
ind4 = where( abs(alog10(tt) - alog10(0.060)) LE 0.007 AND abs(rr-2.05) LE 0.03, ni4)
print,ni4

print,num[ind4],match[ind4],tt[ind4],rr[ind4]



END







